## Import data from .RData
load("data/cont_snotel.RData")

# Import seasonal snow zone polygons
SSZ <- st_read('data/MODIS_Snow_Zones/SSZ_0cc.shp')
## Reading layer `SSZ_0cc' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\SSZ_0cc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357223 ymin: 1163534 xmax: -437223.4 ymax: 3151034
## Projected CRS: NAD_1983_Albers
PSZ <- st_read('data/MODIS_Snow_Zones/PSZ_0cc.shp')
## Reading layer `PSZ_0cc' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\PSZ_0cc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357223 ymin: 1459534 xmax: -437223.4 ymax: 3149702
## Projected CRS: NAD_1983_Albers
# Group sites by site_id
cont_snotel_site <- cont_snow_data %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise()
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# transform files to CRS:4326
snotel_sites <- st_as_sf(cont_snotel_site, coords = c('longitude', 'latitude'), crs = 4326) #fix before next run
PSZ_4326 <- st_transform(PSZ, crs = 4326)

# clip SNOTEL Sites to PSZ
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
PSZ_snotel_sites = st_intersection(PSZ_4326, snotel_sites)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
mapview(PSZ_snotel_sites) + 
  mapview(PSZ_4326)
# Filter all SNOTEL to PSZ SNOTEL Sites Only
PSZ_cont_snotel <- cont_snow_data[cont_snow_data$site_id %in% PSZ_snotel_sites$site_id,]

# PSZ SNOTEL Site, Calculate spring days with new SWE and new SWE depth 
PSZ_cont_snotel_spring <- PSZ_cont_snotel %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))

# SNOTEL Analysis, calculate the number and amount of positive SWE between 
cont_snotel_spring <- cont_snow_data %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))

# PSZ determine average number of days with increased SWE per spring at each site
PSZ_cont_snotel_site <- PSZ_cont_snotel_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years, 
            newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Determine average number of days with increased SWE per spring at each site
cont_snotel_site <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years, 
            newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# View data
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
        layer.name = "AVG Sping Days (All Sites)" ,crs = 4326, grid = FALSE) + 
  mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
        layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE)
# Days of new snow at SNOTEL Sites in Western US
#Summarize the number of new SWE days and SWE depth at all sites per year
cont_snotel_site_yr <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a box plot for all sites
spring_snow_days <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_days)
#Summarize the number of new SWE days and SWE depth at PSZ sites per year
PSZ_cont_snotel_site_yr <- PSZ_cont_snotel_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a boxplot for all sites
PSZ_spring_snow_days <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_days)
summary(cont_snotel_site_yr$elev)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     128    1817    2393    2290    2772    3544
hist(cont_snotel_site_yr$elev)

cont_snotel_site_yr <- cont_snotel_site_yr %>% 
  mutate(elev_band = ifelse(elev >= 2000, "High (>=2500)", 
                            ifelse(elev < 2000 & elev >= 1200, "Mid (1500-2500)", "Low (<= 1500)")))
# Elevation vs new snow count
elev_newsnow <- plot_ly(cont_snotel_site, x = ~avg_spring_days, y = ~elev, color = ~state, type = 'scatter', mode = 'markers', 
                        hoverinfo = 'text',
                        text = ~paste('</br> avg_new_days: ', avg_spring_days,
                                      '</br> elev: ', elev,
                                      '</br> site_name: ', site_name,
                                      '</br> state: ', state)) %>% 
  layout(title = "Average Increased SWE Days by Elevation")
 
elev_newsnow
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
elev_newsnow_depth <- plot_ly(cont_snotel_site, x = ~avg_spring_snow, y = ~elev, color = ~state, type = 'scatter', mode = 'markers', 
                        hoverinfo = 'text',
                        text = ~paste('</br> avg_new_SWE (mm): ', avg_spring_snow,
                                      '</br> elev: ', elev,
                                      '</br> site_name: ', site_name,
                                      '</br> state: ', state)) %>% 
  layout(title = "Average Increased SWE Depth (mm) by Elevation")
 
elev_newsnow_depth
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
elev_boxplot <- ggplot(cont_snotel_site_yr, aes(x = elev_band, y= newsnow_days, color = elev_band)) +
  geom_boxplot()+
  ggtitle("New Snow Events by Elevation")
ggplotly(elev_boxplot)
# Spring Snow depth
spring_snow_depth <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
  geom_boxplot()+
  ggtitle("Avg Depth of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_depth)
# Spring Snow depth
PSZ_spring_snow_depth <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
  geom_boxplot()+
  ggtitle("Depth of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_depth)
# spring_snow_days_yr <- ggplot(cont_snotel_site_yr, aes(x = year, y = newsnow_days, group = state, color = state)) + 
#   geom_line()
# 
# ggplotly(spring_snow_days_yr)